Antibody–nanobody combination increases their neutralizing activity against SARS-CoV-2 and nanobody H11-H4 is effective against Alpha, Kappa and Delta variants

The global spread of COVID-19 is devastating health systems and economies worldwide. While the use of vaccines has yielded encouraging results, the emergence of new variants of SARS-CoV-2 shows that combating COVID-19 remains a big challenge. One of the most promising treatments is the use of not only antibodies, but also nanobodies. Recent experimental studies revealed that the combination of antibody and nanobody can significantly improve their neutralizing ability through binding to the SARS-CoV-2 spike protein, but the molecular mechanisms underlying this observation remain largely unknown. In this work, we investigated the binding affinity of the CR3022 antibody and H11-H4 nanobody to the SARS-CoV-2 receptor binding domain (RBD) using molecular modeling. Both all-atom steered molecular dynamics simulations and coarse-grained umbrella sampling showed that, consistent with the experiment, CR3022 associates with RBD more strongly than H11-H4. We predict that the combination of CR3022 and H11-H4 considerably increases their binding affinity to the spike protein. The electrostatic interaction was found to control the association strength of CR3022, but the van der Waals interaction dominates in the case of H11-H4. However, our study for a larger set of nanobodies and antibodies showed that the relative role of these interactions depends on the specific complex. Importantly, we showed Beta, Gamma, Lambda, and Mu variants reduce the H11-H4 activity while Alpha, Kappa and Delta variants increase its neutralizing ability, which is in line with experiment reporting that the nanobody elicited from the llama is very promising for fighting against the Delta variant.

. (A) Schematic description of the S protein of SARS-CoV-2, which consists of the S1 and S2 subunits. (B) SARS-CoV-2 S protein binds to human ACE2 before entering cells. (C) H11-H4 and CR3022 bind to S protein, preventing the virus from entering cells. The 3D structures of H11-H4 and CR3022 bound to RBD are shown in all-atom (D) and coarse-grained (E) models. Table 1. Dissociation constant K d (nM) obtained by in vitro experiment. The experimental binding free energy ∆G exp was converted from K d using G exp = RTlnK d . Binding free energy ∆G bind (kcal/mol) was obtained using coarse-grained umbrella sampling and Eq. (4) for the H11-H4-RBD, CR3022-RBD, and H11-H4 + CR3022-RBD complexes. Shown is the WT case . www.nature.com/scientificreports/ which is greater than K d obtained by Tian et al. 20 for CR3022, suggesting that H11-H4 binds to RBD weaker than the CR3022 antibody. However, when comparing with K d reported by Yuan et al. 21 (Table 1), we see that H11-H4 binds to RBD more strongly than CR3022. To solve this dispute we will calculate binding affinity using molecular simulation. It is important to note that nanobodies can be used alone or in combination with antibodies in the treatment of severely ill patients with Covid-19 ( Fig. 1C) 25 . The binding affinity of antibodies to SARS-CoV-2 was computationally studied 26,27 , but the binding free energy of nanobodies has not been calculated although their interaction with RBD was explored using molecular modeling. Moreover, the effect of the combination of antibodies and nanobodies on their neutralizing ability has not been theoretically investigated. Therefore, in this paper, using the coarse grained model and umbrella sampling, we will calculate the binding free energy of the H11-H4 nanobody with RBD and study how the combination of the CR3022 antibody and the H11-H4 nanobody changes their ability to neutralize SARS-CoV-2.

Material and methods
PDB structures of the three studied systems. The structures of H11-H4-RBD, CR3022-RBD, and H11-H4 + CR3022-RBD complexes were extracted from the Protein Data Bank with PDB ID: 6ZH9 25 . Modeler package 42 was used to add the missing residues. The structure of H11-H4 + CR3022-RBD complex is shown in Fig. 1D (all-atom) and Fig. 1E (coarse-grained) prepared by using the PyMOL package 43 . All mutations including variants Alpha, Beta, Gamma, Kappa, Delta, Lambda and Mu were generated by using the mutagenesis tool in PyMOL package.
All-atom molecular dynamics simulations. All-atom molecular dynamics (MD) simulations were performed using the CHARMM36M force field 44 implemented in the GROMACS 2016 package 45 at 310 K and isotropic pressure of 1 bar, which was obtained using the v-rescale 46 and Parrinello-Rahman 47 algorithms, respectively. The water model TIP3P 48 was used for all systems. Bond lengths were constrained by the linear constraint solver (LINCS) algorithm 49 , allowing a time step of 2 fs.
Electrostatic and van der Waals interactions were calculated with a cutoff of 1.4 nm, and the non-bonded interaction pair-list was updated every 10 fs. The Particle Mesh Ewald algorithm 50 was used to treat long-range electrostatic interactions. Periodic boundary conditions were applied in all directions. The energy of the system was first minimized by using the steepest-descent algorithm, then a short 3 ns MD simulation was performed in the NVT and NPT ensembles. Production MD simulation of 100 ns was performed with the help of the leapfrog algorithm 51 . For each complex, using the "gmx_mpi cluster" tool available in GROMACS, we grouped the snapshots collected from the 100 ns of conventional MD simulation into clusters. We then selected 5 representative structures from the five most populated clusters and used them as the initial configuration for running 5 trajectories of steered molecular dynamics (SMD) simulations [52][53][54][55] .
Steered molecular dynamics. We carried out SMD simulations to pull H11-H4 or CR3022 from the binding region of RBD as well as pulling RBD from the binding region of H11-H4 and CR3022 (Fig. 2). In the case of H11-H4-RBD and CR3022-RBD, an external force is applied to a dummy atom, which is linked to the Cα atom closest to the center of mass (COM) of H11-H4 or CR3022. The pulling direction is parallel to the vector connecting COMs of RBD and nanobody or antibody ( Fig. 2A,B). In order to prevent RBD from drifting under the action of an external force, its backbone was restrained, but the side chain could fluctuate. The choice of pulling direction is different in the case of H11-H4 + CR3022-RBD due to the presence of three molecules. In this case an external force is applied to a dummy atom that is bonded to the Cα atom closest to the COM of RBD, and the pulling direction is along the line connecting The RBD COM in perpendicular to the line connecting the COMs of H11-H4 and CR3022 (Fig. 2C). During the SMD simulation the backbone of H11-H4 and CR3022 was restrained. For convenience, three complexes H11-H4-RBD, CR3022-RBD and H11-H4 + CR3022-RBD were rotated so that the pulling direction was always along the z-axis (Fig. 2).
One of the limitations of unidirectional pulling is that not all rotational states of proteins can be sampled. However, as shown in previous works 26,55 , this approach provides reasonable results on the relative binding affinity of protein-protein complexes.
The pulling force experienced by a stretched molecule is calculated as follows: where k is the stiffness of the spring, v is the pulling velocity, z is the displacement of a real atom connected to the spring in the direction of pulling, respectively. The spring constant k was set to 600 kJ/(mol nm 2 ) (≈ 1020 pN/nm), which is a typical value used in atomic force microscopy (AFM) experiments 56 . www.nature.com/scientificreports/ Using the force-displacement profile obtained from SMD simulations, the non-equilibrium work (W) performed by the pulled chain (H11-H4, CR3022 or H11-H4 + CR3022) was estimated using the trapezoidal rule: where N is the number of simulation steps, F i and z i are the force determined by Eq. (1) and the position at step i, respectively.
To estimate the binding free energy (∆G), we can use Jarzynski's equality 57 extended to the case when the applied external force grows at a constant speed v 58 : here . . . N is the average over N trajectories, z t is the time-dependent displacement, and W t is the non-equilibrium work at time t determined by Eq. (2). From Eq. (3), we can extract the equilibrium free energy if the number of SMD trajectories is large enough and the pulling is sufficiently slow. Therefore, this approach is practical for small systems 59 but not for large systems such as those studied in this work. However, we can estimate the non-equilibrium binding and unbinding barriers separating the transition state (TS) from the bound state at t = 0 and the unbound state at t end 60 , which allows us to discern weak binding from strong binding.
Rectangular boxes with dimension of 10 × 9 × 25 nm 3 , 10 × 11 × 25 nm 3 and 10 × 18 × 25 nm 3 were used for H11-H4-RBD, CR3022-RBD and H11-H4 + CR3022-RBD, respectively. The complexes were immersed in a 0.15 M sodium chloride salt solution and counter ions were added to neutralize the system. In order to show that our results on the relative binding affinity is independent of the pulling speed, for each system, 5 different trajectories were run at v = 0.5 and 1 nm/ns. Coarse-grained simulation. Since the combination of SMD with Jarzynski's equality does not allow us to calculate the equilibrium binding free energy, we will use umbrella sampling (US). However, this approach is very time consuming if we use all-atom models, because in our case the proteins and antibodies are large. Therefore, we used the MARTINI 2.2 force field developed for coarse-grained (CG) modeling of biological systems such as biological membrane, protein, nucleotide, and etc [61][62][63] . This force field is accurate enough for extracting the interaction energy for a pair of proteins in an aqueous environment from constraint force profiles. The standard MARTINI water model was used with a minimum distance between water beads of 1.0 nm 64 . The system was neutralized by adding sodium chloride salt solution. The temperature was set at T = 300 K with a Berendsen thermostat, and pressure was set at p = 1.0 bar with a Berendsen barostat 65,66 . Bond lengths in the aromatic amino acid side chains and the bonds between the backbone and side chains were constrained by the LINCS algorithm 49 . www.nature.com/scientificreports/ To perform coarse-grained umbrella sampling 67 (CG-US) simulations, we made a series of configurations along the z-axis involving 81 windows each of 0.1 nm (Fig. S1). Here z is the reaction coordinate (RC). The choice of the z-axis has been already described in the SMD method. Namely, for CR3022-RBD and H11-H4-RBD this axis connects two COMs (Fig. S1A,B), while for H11-H4 + CR3022-RBD it is parallel to the line connecting COMs of H11-H and CR3022 (Fig. S1C).
To create an initial configuration for the first window, energy minimization was performed and the neutralized and solvated structure was simulated for 1 ns with position restraints throughout the structure to allow the solvent to equilibrate around the solute. Temperature and pressure were relaxed for 10 ns. The resulting conformation was then used as the initial conformation in a subsequent 100 ns run without position restraints. The last snapshot obtained in this run will be used as the initial configuration for the first window in CG-US simulations.
To generate the initial configuration for other windows, we pulled antibody, nanobody or RBD to the corresponding window. Then we performed energy minimization and equilibration using a 5 ns MD simulation restraining the distance between COMs of subsystems. The last snapshot obtained in this simulation will be used as an initial conformation for the production run.
To hold one chain (H11-H4, CR3022 or RBD) around the center of each window, we applied a bias harmonic potential with a spring constant of 600 kJ/mol/nm 2 to make sure that the interacting surface of both targets is not change. To get a good sampling, for each window, we performed a conventional MD production run of 1000 ns. The WHAM procedure 68 is then used to determine a one-dimensional potential of mean force (1D PMF) as a function of the reaction coordinate z.
The binding free energy ( G bind ) is defined as the difference between the free energies in the bound and unbound states 69 : here G 1D (z) is the 1D PMF as a function of z, k B is the Boltzmann constant, and T is the absolute temperature. Symbols ∫ bound and ∫ unbound refer to summation over bound and unbound regions, respectively. To determine the cut-off distance between the bound and unbound states we calculated the number interchain contacts as a function of the distance between pulled and nonpulled chains in CG-US simulations. Then the cutoff distance is the distance above which interchain contacts disappear (Fig. S2).

Measures used in data analysis. A hydrogen bond (HB) is formed if the distance between donor D and
acceptor A is less than 0.35 nm, the H-A distance is less than 0.27 nm, and the D-H-A angle is greater than 135 degrees. A non-bonded contact (NBC) between two residues is formed if the shortest distance between their atoms is within 0.39 nm. 2D contact networks of HBs and NBCs of CR3022-RBD and H11-H4-RBD were displayed using the LIGPLOT package 70 . The standard deviation (Er) are approximately expressed as follows: where N is the total number of data points in the data set, x i is the individual value of the ith in the data set, and x is the mean value of the data set.

Results and discussion
Hydrogen bond and non-bonded contact networks of CR3022-RBD and H11-H4-RBD complexes: analysis based on the PDB structure. Using the 6ZH9 PDB structure, we build networks of hydrogen bonds (HBs) and non-bonded contacts (NBCs) of H11-H4 and CR3022 with RBD ( Fig. S3A-D). The numbers of H11-H4 and CR3022 residues that form HB and NBC with RBD are 11 and 19, respectively. There are 9 and 10 HBs for H11-H4-RBD and CR3022-RBD, respectively, while the numbers of NBCs of H11-H4-RBD and CR3022-RBD correspond to 14 and 20. The number of HBs and NBCs in the crystal structure cannot determine the binding affinity, since other factors also matter. However, more HBs and NBCs may indicate higher binding affinity, which suggests that CR3022 has a higher binding affinity for RBD than H11-H4. To verify this we will carry out SMD and coarse-grained umbrella simulations.
Binding affinity of H11-H4 and CR3022 to RBD: SMD results. CR3022 binds to RBD more strongly than H11-H4 and combination of antibody and nanobody enhances their neutralizing activity. Force-time profiles obtained at v = 0.5 nm/ns for the three complexes (Fig. 3A, Table 2) show that CR3022 (F max = 1214.2 ± 21.2 kcal/ mol) binds to RBD more strongly than H11-H4 (F max = 925.6 ± 15.2 kcal/mol) to RBD. It should be noted that the rupture force F max appears to be quite high due to the fast pulling. In the so-called Bell approximation, where the transition state separating the bound state from the unbound state is independent of external force, F max ~ ln(v) 71 , where v is the puling speed. Beyond the Bell approximation, the dependence of F max on v is more complex 72 .
As expected, F max increases with increasing of pulling speed (Tables 2 and S1, Figs. 3A and S4A). The unbinding time t max of CR3022-RBD is also longer than H11-H4-RBD, and this time decreases with increasing v. It is important to note that if RBD is simultaneously extracted from the CR3022 antibody and H11-H4 nanobody, then at v = 0.5 nm/ns, F max = 2034.9 ± 27.7 kcal/mol is required, which is approximately twice as much as in CR3022-RBD and H11-H4-RBD. Therefore, the combination of nanobody and antibody is expected to increase the binding affinity for RBD, which increases their neutralizing activity. www.nature.com/scientificreports/ Since the non-equilibrium work (W) is better than F max for characterizing the relative binding affinity 73 , we will look at it in detail. W increased rapidly until the pulled molecule (H11-H4, CR3022 or RBD) left the binding region, reaching a stable value when two subsystems ceased to interact (Fig. 3B), and also increases with v 74 . For v = 0.5 nm/ns, we obtained W = 101.6 ± 3.4, 208.6 ± 5.3 and 461.3 ± 5.7 kcal/mol for H11-H4-RBD, CR3022-RBD, and H11-H4 + CR3022-RBD, respectively (Table 2). Therefore, similar to F max , our results obtained for W further support the fact that CR3022 is more active than H11-H4 and their combination increases their binding strength to RBD. Figure 3C displays the time dependence of the non-equilibrium binding free energy (∆G) estimated from Eq. (3) for the three complexes at v = 0.5 nm/ns. The maximum corresponds to the transition state (TS) with ∆G = ∆G TS . We have ∆G bound = ∆G(t 0 = 0) ≈ 0 kcal/mol at the beginning of the bound state, while the unbound state occurs at the end of simulation ∆G unbound = ∆G(t end ) ≈ 0 kcal/mol. Thus, the binding and unbinding free energy energies (barriers), defined as ∆∆G bind = ∆G TS − ∆G unbound and ∆∆G unbind = ∆G TS − ∆G bound , are roughly equal. ∆∆G unbind = 82.5 ± 2.1, 140.2 ± 2.9, and 379.0 ± 4.2 kcal/mol for H11-H4-RBD, CR3022-RBD, and H11-H4 + CR3022-RBD, respectively, and ∆∆G bind = 82.1 ± 2.5, 137.9 ± 3.7 and 378.5 ± 4.4 kcal/mol ( Table 2). This provides further evidence that CR3022 binds to RBD more tightly than H11-H4, and the binding affinity is higher if both H11-H4 and CR3022 are combined.
To ensure that our result does not depend on the pulling speed, we also conducted SMD simulations for v = 1 nm/ns. Although F max , W, and the non-equilibrium binding free energy increase with increasing v, the main conclusion about the relative binding affinities of the three complexes remains the same (Fig. S4, Table S1).
Therefore, our SMD data indicate that CR3022 binds more strongly to RBD than H11-H4, which is consistent with the experiment of Tian et al. 20 and Huo et al. 25 (Table 1). Measuring K d , Yuan et al. 21 reported that the binding affinity of CR3022 for RBD is lower than that reported by Huo et al. 25 for H11-H4. From this point of view our  Table 2. Rupture force (F max ), rupture time (t max ), work of the external force (W), and non-equilibrium binding free energy ∆G bind and ∆G unbind for the H11-H4-RBD, CR3022-RBD, and H11-H4 + CR3022-RBD complexes. The results were obtained from five independent SMD trajectories for the WT case at a pulling speed of v = 0.5 nm/ns. The errors represent standard deviations. www.nature.com/scientificreports/ result is in conflict with Yuan et al. 21 and Huo et al. 25 The discrepancy may be caused by different techniques used by the two groups. Namely, Yuan et al. 21 used biolayer interferometry binding assays, while isothermal titration calorimetry was employed by Huo et al. 25 The advantage of our computational study is that we used the same model to compare the relative binding affinity, giving us confidence that CR3022 is a better binder than H11-H4.

Binding of H11-H4 to RBD is driven by vdW interaction, but binding of CR3022 and H11-H4 + CR3022 is driven by electrostatic interactions.
A cutoff of 1.0 and 1.2 nm for van der Waals and electrostatic energies was applied to investigate the interaction of H11-H4-RBD, CR3022-RBD and H11-H4 + CR3022-RBD complexes. Fig. 4A1,A2 display the time dependence of the total non-bonded interaction energy E total , which is the sum of electrostatic (E elec ) and van der Waals (E vdW ) energies of H11-H4, CR3022, and H11-H4 + CR3022 interacting with RBD. These results were averaged over five SMD trajectories.
At v = 0.5 nm/ns, when bound, the E elec of H11-H4-RBD, CR3022-RBD and H11-H4 + CR3022-RBD starts at a negative value, but in all three cases E elec eventually reaches ≈ 50 kcal/mol in the unbound state. E vdW of three complexes is also negative in the bound state reaching 0 kcal/mol in the unbound state (Fig. 4A2). Neglecting the contribution of entropy, the results shown in Fig. 4A1,A2 reaffirm the ordering of stability H11-H4 + CR3022 > CR3022 > H11-H4.
We calculated the mean interaction energy in the bound state by averaging over the time window [0, t max ], where t max is shown in Table 3. At v = 0.5 nm/ns, for CR3022-RDB, we obtained E elec = − 252.8 ± 3.7 kcal/mol, which is clearly lower than E vdW = − 77.1 ± 1.3 kcal/mol, implying that binding of CR3022 to RBD is driven by electrostatic interactions. This observation was also obtained previously 26 . The opposite occurs for the case of H11-H4, where the vdW interaction (E vdW = − 61.8 ± 1.2 kcal/mol) is lower than the electrostatic interaction (E elec = − 8.9 ± 0.7 kcal/mol), indicating that the vdW interaction dominates, but not the electrostatic interaction. Thus, the nature of binding of the H11-H4 nanobody is very different from CR3022 and the question of whether this remains true for other nanobodies is left for future research.
If H11-H4 and CR3022 simultaneously bind to RBD, we obtained E elec = − 355.5 ± 2.3 kcal/mol and E vdW = − 146.1 ± 1.1 kcal/mol, which means that as in the single CR3022 case, the electrostatic interaction is more important than the vdW interaction in stabilizing the complex with RBD. The role of electrostatic and vdW interactions revealed in SMD simulations with v = 0.5 nm/ns remains unchanged for other pulling speed (v = 1 nm/ns) (Fig. S5, Table S2).

Role of specific residues in binding of H11-H4 and CR3022 to RBD.
To understand the role of each residue at the interface (Fig. 4B1,B2) in stabilization of the three complexes, we calculated its interaction energy in the [0, t max ] time window for pulling speed v = 0.5 nm/ns. For CR3022-RBD, residues Lys378(C), Lys386(C) and Asp428(C) of RBD, and residues Asp56(A) and Glu58(A) of CR3022 have the total non-bonded interactions smaller than − 20 kcal/mol (Fig. 4B1). For H11-H4-RBD, residues Glu484(C) and Gln493(C) of RBD and residue Arg52(D) of H11-H4 have the interaction energy smaller than − 20 kcal/mol (Fig. 4B2) The role of electrostatic and vdW interactions in the binding of nanobodies and antibodies to RBD depends on the specific system. Our previous work 26 showed that the electrostatic interaction governs the binding of CR3022 to RBD, while in the present work the vdW interaction is found to be more important for H11-H4 nanobody. An interesting question emerges is if this conclusion is valid for other systems. To answer this question, we calculated the interaction energy for the ten antibody-RBD complexes and ten nanobody-RBD complexes using their PDB structures and the CHARMM36M force field with the TIP3P water model.
For antibodies, electrostatic interaction dominates over vdW interaction for the five antibodies, while vdW interaction takes over electrostatic interaction for the other five antibodies (Table S3). For nanobodies, the vdW interaction is more important than the Coulomb interaction in five cases, while the opposite occurs in the other four complexes. In the case of WNb 10-RBD, their role is almost the same (Table S4). Consequently, which interaction is dominant in the association of the antibodies and nanobodies with the spike protein depends on the specific system. Effects of mutations on binding affinity of H11-H4 to RBD: SMD results. As mentioned in the previous section, for the WT case, residue 484 makes an important contribution to the stability of the H11-H4-RBD complex. It has recently been demonstrated that this residue decreases the neutralizing activity of antibodies and nanobodies against the Covid-19 variants 37,41 (see Table S5 for mutation points in some variants).  (Table S5). Note that CR3022 does not have contact with all R346, K417, L452, T478, E484, F490 and N501 residues of RBD, where the mutation is made for the aforementioned SARS-CoV-2 variants (Fig. S6). Therefore, we carried out SMD simulation only for H11-H4-RBD.
Beta, Gamma, Lambda and Mu variants reduce the binding affinity of H11-H4 to RBD. As seen from Fig. 5A1 Table 5). The strong attractive interaction becomes even repulsive after E484K mutation. Meanwhile, the decrease in the interaction of H11-H4 with the Lambda variant occurs predominantly due to L452Q and F490S, which change the total non-bonded interaction energy from − 0.2 to − 8.4 kcal/mol (WT)  Alpha, Kappa and Delta variants increase the binding affinity of H11-H4 to RBD. In addition to the Beta, Gamma, Lambda and Mu variants, we also examined the binding affinity of H11-H4 to Alpha, Kappa and Delta Table 3. Non-bonded interaction energies (kcal/mol) of the H11-H4-RBD, CR3022-RBD, and H11-H4 + CR3022-RBD complexes. The results were obtained for a [0-t max ] time window and averaged from five SMD trajectories performed at a pulling speed of v = 0.5 nm/ns. The errors represent standard deviations.   Table 4), implying that H11-H4 can neutralize these variants better than WT. For the Alpha variant, although the mutation point N501Y does not significantly contribute to the stability of H11-H4-RBD (Table 5), the binding affinity is insignificantly stronger than that of WT (Fig. 5B1-B3; Table 4). The total non-bonded interaction energy of N501Y slightly drops from 0 (WT) to − 0.1 kcal/mol. For the Kappa variant, the E484Q mutation destabilizes the H11-H4-RBD complex, as the corresponding total non-bonded interaction energy increases from − 156.3 (WT) to − 20.1 kcal/mol ( Table 5). The L452R mutation also weakens the interaction with H11-H4 due to an increase in total non-bonded interaction energy from − 0.2 (WT) to 58.6 kcal/mol (Table 5). Based on the total non-bonded interaction energy obtained at mutation positions 484 and 452, we cannot explain why the Kappa variant enhances the stability of H11-H4-RBD complex. Same as the Kappa variant, the L452R mutation of the Delta variant has an increase in the total non-bonded interaction energy from − 0.2 (WT) to 56.5 kcal/mol (Table 5), but the binding affinity is still much higher than WT. So what is the reason for the increased binding affinity between H11-H4 and RBD in the Kappa and Delta variants?
To solve these issues, we calculated the total interaction energy not only for the residues related to the mutation points, but also for all important residues (Fig. 6A,B). For WT, the total energy is − 89.3 kcal/mol, which is higher than Alpha (− 99.1 kcal/mol), Kappa (− 118.4 kcal/mol) and Delta (− 129.3 kcal/mol). For Gamma, Mu, Beta and Lambda we obtained 368.2, 335.1, 320.9 and − 49.3 kcal/mol, respectively, which is clearly higher than for WT. Therefore, the order of stability is as follows: Gamma < Mu < Beta < Lambda < WT < Alpha < Kap pa < Delta. This finding is consistent with a report that nanobodies elicited from a llama could neutralize the Delta variant 41 . In addition, we also predict that H11-H4 maybe an excellent candidate to treat Alpha and Kappa variants.
Binding free energy of H11-H4 and CR3022 to RBD: coarse-grained umbrella sampling results. The SMD, known as a method used to investigate the unbinding process of a small molecule to other molecules, is capable of predicting relative binding affinity but cannot be used to calculate the binding free energy. Overall, although the SMD method has provided a good correlation with experimental results 26,52,55,75 , their predictions are not always perfect. Therefore, we also used coarse-grained umbrella sampling to determine the binding free energies in an effort to elucidate the interactions of H11-H4-RBD, CR3022-RBD and H11-H4 + CR3022-RBD complexes.
The MARTINI CG-US was used to estimate the binding free energy (∆G bind ). To show that the equilibrium phase has been reached, we calculated the 1D PMF for three time intervals of 500, 800 and 1000 ns. Since the Table 4. Rupture force (F max ), pulling work (W), and non-equilibrium binding (∆G bind ) and unbinding (∆G unbind ) free energies obtained from five independent SMD trajectories with v = 0.5 nm/ns for H11-H4-RBD. Results are shown for WT and variants Alpha, Beta, Gamma, Kappa, Delta, Lambda and Mu. The errors represent standard deviations.   www.nature.com/scientificreports/ 1D PMF profiles for these windows are essentially the same (Fig. S7) our data was equilibrated. Therefore, the profile obtained from the largest window ( Fig. 7) was used for estimating the binding free energy. In order to use Eq. (4) to extract the binding free energy from the 1D PMF profiles we must estimate the cutoff distance between bound and unbound states, which was defined in "Material and methods". From the distance dependence of the number of interchain contacts (Fig. S2) we obtained the cutoff distance of 4.9, 2.8 and 3.0 nm for H11-H4 + CR3022-RBD, H11-H4-RBD and CR3022-RBD, respectively. Using these cutoff distances, Eq. (4) and the 1D PMF profiles shown in Fig. 7, we obtained ∆G bind = − 19.8 kcal/mol for H11-H4-RBD and − 21.4 kcal/ mol for CR3022-RBD ( Table 1). The low value of ∆G bind of CR3022 indicates that this antibody tightly binds to the S protein which is consistent with previous computational studies 26 . Moreover, in agreement with the SMD results shown above and the in vitro results reported by Huo et al. 25 (K d = 11.8 nM, ∆G bind = − 10.9 kcal/mol for H11-H4) and Tian et al. 20 (K d = 6.3 nM, ∆G bind = − 11.3 kcal/mol for CR3022), H11-H4 binds to RBD weaker than CR3022. Obviously, the ∆G bind value obtained with the MARTINI CG-US is much lower than the experimental data, which may be related to the force field we used and the complexity of the studied systems.
According to our CG-US results, the relative binding affinity of CR3022 and H11-H4 is ∆G bind (H11-H4)/∆G bind (CR3022) = − 19.8/− 21.4 = 0.93, which is not too far from experimental value ∆G exp (H11-H4)/∆G exp (CR3022) = − 10.9/− 11.3 = 0.97. Thus, although the difference between simulation and experiment in absolute binding free energy is quite large, agreement on the relative binding affinity is satisfactory. In addition, bearing in mind that it is very challenging to obtain the absolute binding free energy even for small ligands using different MD based methods, our results are reasonable.
For H11-H4 + CR3022-RBD we obtained ∆G bind = − 23.9 kcal/mol, which shows that, consistent with the SMD results, the combination of antibody and nanobody enhances their binding affinity. Since

Discussion and conclusion
Combining various computational methods, we studied the association of H11-H4, CR3022, and both H11-H4 and CR3022 with RBD. A number of interesting results have been obtained. SMD simulation showed that the H11-H4 nanobody binds to RBD weaker than CR3022, which is consistent with the binding free energy ∆G bind computed from the coarse-grained US. Our theoretical estimates of the binding affinity are in good agreement with the experimental results presented by Tian et al. 20 and Huo et al. 25 for H11-H4 and CR3022 interacting with SARS-CoV-2, but for the CR3022-RBD complex they contradict Yuan et al. 21 (Table 1). To clarify the contradiction, recall that Tian et al. 20 reported a K d of 6.3 nM for CR3022-RBD, which is lower than Kd = 115 nM of Yuan et al. 21 (Table 1). The difference in binding affinity may be due to differences in experimental conditions of these two groups, as discussed by Yuan et al. 21 . Namely, CR3022 was expressed as an Fab in Yuan et al. 21 , but singlechain fragment variable (scFv) in Tian et al. 20 ; CR3022 was expressed in mammalian cells in Yuan et al. 21 , but in E. coli in Tian et al. 20 . SARS-CoV-2 RBD was expressed in insect 150 cells in Yuan et al. 21 , but in mammalian cells in Tian et al. 20 . Since the result of Yuan et al. 21 does not agree with our simulations, we assume that their experimental conditions do not match our modeling. However, we cannot confirm this with molecular simulations because none of the existing models is able to capture differences in the behavior of proteins in different cells. We predict that the concurrent binding of H11-H4 and CR3022 to RBD results in a higher binding affinity than when they are individually associated with RBD. Thus, the combination of H11-H4 and CR3022 enhances the neutralization of SARS-CoV-2, and this could open up a new treatment strategy for Covid-19. Whether this conclusion holds for the other antibody-nanobody pairs is a matter of further clarification.
Stability of the H11-H4-RBD complex is mainly contributed by the vdW interaction, while electrostatic interaction is more important for the CR3022-RBD and H11-H4 + CR3022-RBD complexes. In general, the role of vdW and electrostatic interaction in the binding of antibodies and nanobodies to SARS-CoV-2 depends on the specific case.
Our computational study found CR3022 to be a better candidate for treating Covid-19 than H11-H4, but only for WT. It is important to note that H11-H4 shows a high ability to neutralize the Alpha, Kappa and the highly dangerous Delta variants, and this fact is consistent with recent experiment.
In our work, the advantage of SMD is that it can provide all-atom description, but the disadvantage is that it does not allow the calculation of the free energy at equilibrium. In contrast, the equilibrium binding free energy can be obtained by coarse-grained umbrella sampling, but only at the coarse level. Nevertheless, these two approach complement each other, leading to a satisfactory description of the experiment.
For protein-ligand systems, SMD has been shown to be as efficient as the MM-PBSA method but computationally faster due to fast pull 76 . Although a similar analysis has not been carried out for protein-protein interactions, in conjunction with previous woks 26,55,75 , the present study shows that this method is useful for characterizing the relative binding affinity of protein-protein complexes. It would be useful to compare the SMD www.nature.com/scientificreports/ method with the MM-PBSA and other MD-based methods for estimating the absolute binding free energy of these complexes. In the reaction coordinate space, PMF shown in Fig. 7 is one-dimensional (1D), and mapping the multidimensional free energy landscape to the 1D profile is an approximation. However, Z in Fig. 7 is the radial distance in real 3D space, which can reflect the 3D nature of the problem. This may be one of the reasons why umbrella sampling is one of the best methods for calculating free energy 77 . In other words, 1D PMF is adequate for our problem.

Data availability
The data files are available in the "SciRep_data.zip" file, which includes: SMD simulation data: They are located in the "SMD_data" folder. The "WT" subfolder contains binding affinity and interaction energy data obtained from SMD simulation for the wild type. The "MUTATION" subfolder contains data on Covid-19 variants. Coarsegrained simulation data: The data obtained from umbrella sampling coarse-grained simulations are presented in the "US-CG" folder.